#include "../bulk_storage/bulk_storage.h"
#include "atom.h"
#include "residue.h"
//#include "../aa_neighborhood/aa_neighborhood.h"
#include "aa_table.h"
#include "atom_weight.h"
#include "../metric/xyz.h"

#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>

#include <assert.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void mem_res_print(const struct mem_res* r){
  char s[5];
  int i=1;
  s[4] = '\0';
  for(i=1; i<=r->n; i++){
    memcpy(s, r->atoms[i]->name, 4*sizeof(char));
    printf("%d: %s\n", r->atoms[i]->serialNumber, s);
  }
}


void mem_res_add_atom(struct mem_res* r, const struct atom* a){

  const int INIT_NO_ATOMS = 27;

  /* Alloc 27 pointers if this is the first atom  */
  if(r->n == 0){ 
    r->atoms = malloc(INIT_NO_ATOMS*sizeof(struct atom*));

    /* Alloc more space if we are adding more than 27 atoms  */
 } else if(r->n > INIT_NO_ATOMS){
   printf("realloc() in residue_add_atom().\n");
	fflush(stdin);
    r->atoms = realloc(r->atoms, r->n*sizeof(struct atom*));
 }

  r->n++;
  r->atoms[r->n] = malloc(sizeof(struct atom));
  assert(r->atoms[r->n] != 0);
  memcpy(r->atoms[r->n], a, sizeof(struct atom));
}


void mem_res_load(struct mem_res* mr, 
		  const struct residue* r, 
		  const struct store* s){
  int o=0;
  int r_idx = aa_get_table_offset(r->residue_type);
  const char* asan0 = aa_active_sites[r_idx].first_atom;
  const char* asan1 = aa_active_sites[r_idx].second_atom;

  struct atom* asa0=0;
  struct atom* asa1=0;
  struct atom* a=0;

  mr->n = 0;

  for(o=r->startA; o<=r->stopA; o++){
    a = store_read(s, o);
    mem_res_add_atom(mr, a);
    if( strncmp(asan0, a->name, 4) == 0 ){
      asa0 = a;
    }
    if( strncmp(asan1, a->name, 4) == 0 ){
      asa1 = a;
    }
    store_release(s, o);
  }
  
  mr->res=*r;
  mr->fg = gsl_vector_alloc(3);
  mem_res_compute_functional_group(mr, asa0, asa1);
}


double mem_res_min_effect_dist(const struct mem_res* r0, const struct mem_res* r1){
  double* d0;
  double* d1;
  int i;

  gsl_vector* tmp  = gsl_vector_alloc(3);

  d0 = malloc((r0->n)*sizeof(double));
  assert(d0 != 0);
  d1 = malloc((r0->n)*sizeof(double));
  assert(d1 != 0);

  for(i=1; i<=r0->n; i++){
    atom2vd(tmp, r0->atoms[i]);
    gsl_vector_sub(tmp, r0->fg);
    d0[i] = gsl_blas_dnrm2(tmp);
  }
  for(i=1; i<=r1->n; i++){
    atom2vd(tmp, r1->atoms[i]);
    gsl_vector_sub(tmp, r1->fg);
    d1[i] = gsl_blas_dnrm2(tmp);
  }
  qsort(d0, r0->n, sizeof(double), &compdbl);
  qsort(d1, r1->n, sizeof(double), &compdbl);
  
  if(d0[0] <= d1[0]){return(d0[0]);}
  else {return(d1[0]);}

}


int compdbl(const void* d0, const void* d1){
  if(*((double*)(d0)) < *((double*)(d1))){return(-1);}
  else if(*((double*)(d0)) == *((double*)(d1))){return(0);}
  else{return(1);}
}


double residue_min_effect_dist(const struct residue* r0, 
			       const struct residue* r1, 
			       const struct store* s){
  struct mem_res mr0;
  struct mem_res mr1;
  
  mem_res_load(&mr0, r0, s);
  mem_res_load(&mr1, r1, s);

  return(mem_res_min_effect_dist(&mr0, &mr1));
}


void mem_res_compute_functional_group(const struct mem_res* mr,
				      const struct atom* asa0,
				      const struct atom* asa1){

  gsl_vector* asv0 = gsl_vector_alloc(3);
  gsl_vector* asv1 = gsl_vector_alloc(3);

  assert(mr->fg != 0);
  /* If one atom is null, set it equal to the other.
   * The averaging will take care of the rest. */
  if(asa0 == 0){ asa0 = asa1; }
  if(asa1 == 0){ asa1 = asa0; }

  /* No functional group defined. 
   * Use the sidechain mass center.  */
  if(asa0 == 0){ 
    struct atom** a;
    size_t n;
    mem_res_side_chain(mr, &a, &n);
    mass_center(mr->fg, (const struct atom**)a, n);

  } else {

    /* Take the mean of the two points. */
    atom2vd(asv0, asa0);
    atom2vd(asv1, asa1);
    
    gsl_blas_dscal(0.5, asv0);
    gsl_blas_dscal(0.5, asv1);
    gsl_vector_add(asv0, asv1);

    gsl_vector_memcpy(mr->fg, asv0);
  }
}



void atom2vd(gsl_vector* v, const struct atom* a){
  assert(v != 0);
  gsl_vector_set(v, 0, (double)a->x);
  gsl_vector_set(v, 1, (double)a->y);
  gsl_vector_set(v, 2, (double)a->z);
}


void mass_center(gsl_vector* v, const struct atom** ag, size_t n){

  double w=0;
  double W = 0;
  int i;
  gsl_vector* tmp = gsl_vector_alloc(3);

  assert(v != 0);

  gsl_vector_set_all(v, 0.0);
  
  for(i=0; i<n; i++){
    atom2vd(tmp, ag[i]);
    w = find_atom_mass(ag[i]->element);
    gsl_blas_daxpy(w, tmp, v);
    W += w;
  }
  gsl_blas_dscal(1/W, v);
}


void mem_res_side_chain(const struct mem_res* mr, struct atom*** ag, size_t* n){
  const char main_chain[4][4] = {" N  ",
				 " CA ",
				 " C  ",
				 " O  "};
  const int noatoms = 28;
  int i, j, k=0;
  int done = 0;
  /* Initial space  */
  (*ag) = malloc(noatoms*sizeof(struct atom*));

  /* Look for atoms not on the main chain.  */
  for(i=1; i<mr->n; i++){
    for(j=0; j<4; j++){
      if(strncmp(mr->atoms[i]->name, main_chain[j], 4) == 0){
	done = 1;
      } 
    }
    if(done == 0){
      
      if(k>=noatoms){
	(*ag) = realloc(*ag, (k+1)*sizeof(struct atom*));
      }
      (*ag)[k] = mr->atoms[i];
      k++;
    }
    done =0;
  }
  if(k<noatoms){
    (*ag) = realloc(*ag, k*sizeof(struct atom*));
  }
  *n = k;
}
